# !/usr/bin/env python
# -*- coding: utf-8 -*-
"""
@Time        : 2021/1/2 18:42
@Author      : Albert Darren
@Contact     : 2563491540@qq.com
@File        : orthogonal_polynomial.py
@Version     : Version 1.0.0
@Description : TODO 自己实现勒让德多项式，第一类和第二类切比雪夫多项式及对应切比雪夫多项式零点，拉盖尔多项式，埃尔米特多项式
@Created By  : PyCharm
"""
from numpy import pi
from sympy import expand, cos, Poly, solve, Eq
from sympy.abc import x
import numpy as np


# 自己原创
def legendre_polynomial(symbol_x, degree):
    """
    实现degree次勒让德(Legendre)多项式
    :param symbol_x:符号变量
    :param degree:多项式阶数
    :return:指定degree阶数的勒让德(Legendre)多项式
    """
    legendre_t_list = np.array([1, symbol_x])
    if degree in (0, 1):
        return legendre_t_list[degree]
    for n in range(1, degree + 1):
        legendre_t_list[0] = (2 * n + 1) / (n + 1) * symbol_x * legendre_t_list[1] - n / (n + 1) * legendre_t_list[0]
        # 第一种方式交换legendre_t_list
        legendre_t_list = legendre_t_list[::-1]
        # 第二种方式交换legendre_t_list
        # legendre_t_list[1], legendre_t_list[0] = legendre_t_list[0], legendre_t_list[1]
        # 第三种方式交换legendre_t_list
        # temp = legendre_t_list[1]
        # legendre_t_list[1]=legendre_t_list[0]
        # legendre_t_list[0] = temp
    return expand(legendre_t_list[0])


# 自己原创
def laguerre_polynomial(symbol_x, degree):
    """
    实现degree次拉盖尔(laguerre)多项式
    :param symbol_x:符号变量
    :param degree:多项式阶数
    :return:指定degree阶数的拉盖尔(laguerre)多项式
    """
    laguerre_l_list = np.array([1, 1 - symbol_x])
    if degree in (0, 1):
        return laguerre_l_list[degree]
    for n in range(1, degree + 1):
        laguerre_l_list[0] = (1 + 2 * n - symbol_x) * laguerre_l_list[1] - n * n * laguerre_l_list[0]
        # 第一种方式交换legendre_t_list
        laguerre_l_list = laguerre_l_list[::-1]
        # 第二种方式交换legendre_t_list
        # laguerre_l_list[1], laguerre_l_list[0] = laguerre_l_list[0], laguerre_l_list[1]
        # 第三种方式交换legendre_t_list
        # temp = laguerre_l_list[1]
        # laguerre_l_list[1]=laguerre_l_list[0]
        # laguerre_l_list[0] = temp
    return expand(laguerre_l_list[0])


# 自己原创
def hermite_polynomial(symbol_x, degree):
    """
    实现degree次埃尔米特(hermite)多项式
    :param symbol_x:符号变量
    :param degree:多项式阶数
    :return:指定degree阶数的埃尔米特(hermite)多项式
    """
    hermite_l_list = np.array([1, 2 * symbol_x])
    if degree in (0, 1):
        return hermite_l_list[degree]
    for n in range(1, degree + 1):
        hermite_l_list[0] = 2 * symbol_x * hermite_l_list[1] - 2 * n * hermite_l_list[0]
        # 第一种方式交换hermite_t_list
        hermite_l_list = hermite_l_list[::-1]
        # 第二种方式交换hermite_t_list
        # hermite_l_list[1], hermite_l_list[0] = hermite_l_list[0], hermite_l_list[1]
        # 第三种方式交换hermite_t_list
        # temp = hermite_l_list[1]
        # hermite_l_list[1]=hermite_l_list[0]
        # hermite_l_list[0] = temp
    return expand(hermite_l_list[0])


# 自己原创
def chebyshev_polynomial(symbol_x, degree, polynomial_kind=1):
    """
    实现degree次切比雪夫(Chebyshev)多项式
    :param polynomial_kind: 切比雪夫多项式类别，第一类或者第二类,默认第一类
    :param symbol_x:符号变量
    :param degree:多项式阶数
    :return:指定degree阶数的切比雪夫(Chebyshev)多项式
    """
    if polynomial_kind not in (1, 2):
        raise Exception("Error,chebyshev_polynomial must be the first or second kind,polynomial_kind=1 or 2.")
    chebyshev_t_list = np.array([1, polynomial_kind * symbol_x])
    if degree in (0, 1):
        return chebyshev_t_list[degree]
    for n in range(degree):
        chebyshev_t_list[0] = 2 * symbol_x * chebyshev_t_list[1] - chebyshev_t_list[0]
        # 第一种方式交换chebyshev_t_list
        chebyshev_t_list = chebyshev_t_list[::-1]
        # 第二种方式交换chebyshev_t_list
        # chebyshev_t_list[1], chebyshev_t_list[0] = chebyshev_t_list[0], chebyshev_t_list[1]
        # 第三种方式交换chebyshev_t_list
        # temp = chebyshev_t_list[1]
        # chebyshev_t_list[1]=chebyshev_t_list[0]
        # chebyshev_t_list[0] = temp
    return expand(chebyshev_t_list[0])


# 自己原创
def chebyshev_null_points(degree, a=-1, b=1, polynomial_kind=1):
    """
    利用区间变换和切比雪夫零点求一般区间[a,b]上的插值节点
    :param polynomial_kind: 切比雪夫多项式类别，第一类或者第二类，默认第一类
    :param degree: 切比雪夫多项式次数
    :param a: 区间左端点
    :param b: 区间右端点
    :return:指定区间上的切比雪夫插值节点
    """
    if polynomial_kind == 1:
        return [((b - a) / 2) * cos(((2 * k + 1) / (2 * degree)) * pi) + (b + a) / 2 for k in range(degree)]
    elif polynomial_kind == 2:
        return [((b - a) / 2) * cos((k * pi) / (degree + 1)) + (b + a) / 2 for k in range(1, degree + 1)]
    else:
        raise Exception("Error,chebyshev_polynomial must be the first or second kind,polynomial_kind=1 or 2.")


if __name__ == '__main__':
    # 切比雪夫多项式测试成功，来源详见李庆扬数值分析第5版P61,公式(2.11)
    """
    chebyshev_interpolation_polynomial = Poly(chebyshev_polynomial(x, 6), x)
    # 获得6次切比雪夫多项式最高次项系数
    print(chebyshev_interpolation_polynomial.coeffs()[0])
    print("6次切比雪夫多项式为:{}".format(chebyshev_interpolation_polynomial))
    # 利用T5(x)的零点和区间变换求4次拉格朗日插值节点测试成功，来源详见李庆扬数值分析第5版P64,e.g.4
    print(chebyshev_null_points(5, 0, 1))
    # 用T3(x)的零点做插值点求二次插值多项式，来源详见李庆扬数值分析第5版P94,e.x.11
    print(chebyshev_null_points(3))
    """
    # 测试0-6次Legendre多项式成功，来源详见李庆扬数值分析第5版P61,公式(2.9)
    """
    # 1
    # x
    # 1.5*x**2 - 0.5
    # 2.5*x**3 - 1.5*x
    # 4.375*x**4 - 3.75*x**2 + 0.375
    # 7.875*x**5 - 8.75*x**3 + 1.875*x
    # 14.4375*x**6 - 19.6875*x**4 + 6.5625*x**2 - 0.3125
    for i in range(7):
        print(legendre_polynomial(x, i))
    """
    # 其他常用正交多项式测试成功，来源详见李庆扬数值分析第5版P66
    # 测试第二类切比雪夫多项式零点成功
    """
    print(chebyshev_null_points(7, polynomial_kind=2))
    print(chebyshev_polynomial(x, 7, 2))
    print(solve(Eq(chebyshev_polynomial(x, 7, 2), 0)))
    """
    # 拉盖尔多项式测试成功
    """
    # 1
    # 1 - x
    # x**2 - 4*x + 2
    # -x**3 + 9*x**2 - 18*x + 6
    # x**4 - 16*x**3 + 72*x**2 - 96*x + 24
    for i in range(5):
        print(laguerre_polynomial(x, i))
    """
    # 埃尔米特多项式测试成功
    """
    # 1
    # 2*x
    # 4*x**2 - 2
    # 8*x**3 - 12*x
    # 16*x**4 - 48*x**2 + 12
    for i in range(5):
        print(hermite_polynomial(x, i))
    """
